
#include "ENRadiationHeatFluxUserObject.h"
#include "FEProblem.h"
#include "libmesh/elem.h"
#include "libmesh/fe.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/plane.h"
#include "MooseMesh.h"
#include <fstream>
#include <time.h>
#include "RayLine.h"
#include "MooseRandom.h"

//#include "MemData.h"
using namespace std;

registerMooseObject("RoshanApp", ENRadiationHeatFluxUserObject);
template<>
InputParameters validParams<ENRadiationHeatFluxUserObject>()
{
	InputParameters params = validParams<SideUserObject>();
	params += validParams<RandomInterface>();
	params.addCoupledVar("temperature", "温度场");
	params.addParam< string >("boundary_groups", "The list of boundary groups");
	params.addParam< int >("max_reflect_count", 10, "最大反射次数");
	params.addParam< int >("devide_side", 1, "devide side elemnet");
	params.addParam< vector<Real> > ("transmissivity", "透射率");
	params.addParam< vector<Real> > ("absorptivity", "吸收率");
	params.addParam< bool >("is_seal", true, "is seal");
	params.addParam< bool >("mpi_or_omp", true, "The style of parallel");
	params.addParam< int >("n_threads", 1, "Number of thread while using openMP");
	params.addParam< int >("BSMC", 0, "Deal with BSMC or not");
	params.addParam< bool >("write_or_read", true,  "Write a file of RD, or read RD from a file");
	params.addParam< string >("filename", "The name of file_RD");
	
	return params;
}

ENRadiationHeatFluxUserObject::ENRadiationHeatFluxUserObject(const InputParameters & parameters) :
	SideUserObject(parameters),
	RandomInterface(parameters, *parameters.get<FEProblem *>("_fe_problem"), parameters.get<THREAD_ID>("_tid"), false),
	_temperature(coupledValue("temperature")),
	_max_reflect_count(getParam<int>("max_reflect_count")),
	_devide_side(getParam<int>("devide_side")),
	_transmissivity(getParam< vector<Real> > ("transmissivity")),
	_absorptivity(getParam< vector<Real> > ("absorptivity")),
	_is_seal(getParam< bool >("is_seal")),
	_mpi_or_omp(getParam< bool >("mpi_or_omp")),
	_n_threads(getParam< int >("n_threads")),
	_BSMC(getParam< int >("BSMC")),
	_write_or_read(getParam< bool > ("write_or_read")),
	_filename(getParam< string >("filename"))
{
}

ENRadiationHeatFluxUserObject::~ENRadiationHeatFluxUserObject()
{
	for(unsigned int i=0; i<_all_element.size(); i++)
	{
		delete _all_element[i]->_elem;
	}
}

void ENRadiationHeatFluxUserObject::initialSetup()
{
	string boundary_groups = getParam<string>("boundary_groups");
	vector<BoundaryName> group;
	vector< vector<BoundaryName> > group_set;
	MooseUtils::tokenize<BoundaryName>(boundary_groups, group, 1, ",");
	group_set.resize(group.size());
	for(unsigned int i = 0; i < group.size(); ++i)
	{
		MooseUtils::tokenize<BoundaryName>(group[i], group_set[i], 1, " ");
	}

	for (unsigned int i = 0; i < group.size(); ++i)
	{
		cout << "boundary_group_" << i << ": " << group[i] << endl;
	}

	cout << endl << "Mesh modifying..." << endl;
	
	vector<BoundaryName> boundary = getParam<std::vector<BoundaryName> >("boundary");
	std::set<BoundaryID> boundary_ids;
	for(vector<BoundaryName>::iterator it = boundary.begin(); it != boundary.end(); ++it)
	{
		BoundaryID id = _mesh.getBoundaryID(*it);
		boundary_ids.insert(id);
	}

	MeshBase & mesh = _mesh.getMesh();
	const BoundaryInfo &bnd_info = mesh.get_boundary_info();

	MeshBase::const_element_iterator el = mesh.active_elements_begin();
	const MeshBase::const_element_iterator end_el = mesh.active_elements_end();
	for ( ; el != end_el ; ++el)
	{
		const Elem *elem = *el;

		for (unsigned int side=0; side < elem->n_sides(); ++side)
		{
			if (elem->neighbor_ptr(side))
				continue;

			const Elem *elem_side = elem->build_side_ptr(side).release();
			vector<boundary_id_type> bnd_id_vec;
			bnd_info.boundary_ids(elem, side, bnd_id_vec);
			int bnd_id = bnd_id_vec[0];
			//int bnd_id = bnd_info.boundary_id(elem, side);

			int bnd_in_which_group =-1;
			for(unsigned int i = 0; i < group.size(); ++i)
			{
				for(unsigned int j = 0; j < group_set[i].size(); ++j)
				{
					if(_mesh.getBoundaryID(group_set[i][j]) == bnd_id)
					{
						bnd_in_which_group = i;
						break;
					}
				}
			}

			if(find(boundary_ids.begin(), boundary_ids.end(), bnd_id) == boundary_ids.end())
			{
				delete elem_side;
				continue;
			}

			unsigned int dim = _mesh.dimension();
			FEType fe_type(Utility::string_to_enum<Order>("CONSTANT"), Utility::string_to_enum<FEFamily>("MONOMIAL"));
			FEBase * _fe_face = (FEBase::build(dim, fe_type)).release();
			QGauss * _qface = new QGauss(dim-1, FIRST);
			_fe_face->attach_quadrature_rule(_qface);
			_fe_face->reinit(elem, side);
			const std::vector<Point> normals = _fe_face->get_normals();
			_all_element.push_back(new ENSideElement(elem_side, normals[0], _devide_side, _transmissivity[bnd_in_which_group], _absorptivity[bnd_in_which_group]));

			delete _fe_face;
			delete _qface;
		}
	}

	for (unsigned int i = 0; i < _all_element.size(); i++)
	{
		_center_point_list.push_back(_all_element[i]->_elem->centroid() );
		_epsilon_list.push_back(_all_element[i]->_absorptivity );
	}

	cout << "Mesh modify completed" << endl;

	cout << endl << "Count of boundary element: " << _all_element.size() << endl << endl;

	_temperature_pow4_bar.resize(_all_element.size(), 0);
	_temperature_pow3_bar.resize(_all_element.size(), 0);
	_temperature_bar.resize(_all_element.size(), 0);
	_flux_radiation.resize(_all_element.size(), 0);
	_flux_radiation_jacobi.resize(_all_element.size(), 0);

	clock_t start_time,end_time;
	start_time = clock();

	if (_write_or_read)
	{
		_communicator.barrier();

		for (unsigned int ii=0;ii<_all_element.size();ii++)
		{
			_all_element[ii]->_local_RD.resize(_all_element.size(), 0);
		}

		if(_mpi_or_omp)
		{
			int local_start_id;
			int local_end_id;
			int all_element_size = _all_element.size();
			int quotient = all_element_size / n_processors();
			int remainder = all_element_size % n_processors();
			if(processor_id() < remainder)
			{
				local_start_id = (quotient + 1) * processor_id();
				local_end_id = (quotient + 1) * (processor_id() + 1);
			}
			else
			{
				local_start_id = quotient * processor_id() + remainder;
				local_end_id = quotient * (processor_id() + 1) + remainder;
			}

			//////////////////////开始计算///////////////////////////
			cout << "Computing radiation view factor start" << endl;

			for (int i = local_start_id; i < local_end_id; ++i)
			{
				Real Ai = _all_element[i]->_elem->volume();
				int i_sep = _all_element[i]->_Disintegrate_centre.size();
				Point i_normal = _all_element[i]->getSideElementNormal();
				_all_element[i]->_local_RD.resize(all_element_size, 0.0);
				for (int j = 0; j < all_element_size; ++j)
				{
					Real VF_ij = 0.0;
					if (i == j)
					{
						VF_ij = 0.0;
					}
					else
					{
						Point j_normal = _all_element[j]->getSideElementNormal();
						int j_sep = _all_element[j]->_Disintegrate_centre.size();
						for (int i_in = 0; i_in < i_sep; ++i_in)
						{
							Point centre_i = _all_element[i]->_Disintegrate_centre[i_in];
							Real aera_i = _all_element[i]->_Disintegrate_aera[i_in];
							for (int j_in = 0; j_in < j_sep; ++j_in)
							{
								Point centre_j = _all_element[j]->_Disintegrate_centre[j_in];
								Real aera_j = _all_element[j]->_Disintegrate_aera[j_in];
								if ((_all_element[j]->getSideElementNormal()) * ((centre_j - centre_i).unit()) > (-RaylineZero))
								{
									continue;
								}
								else if ((_all_element[i]->getSideElementNormal()) * ((centre_i - centre_j).unit()) > (-RaylineZero))
								{
									continue;
								}
								else if (Is_Shelter_Engineer(i, j, centre_i, centre_j))
								{
									continue;
								}
								else
								{
									Real cosi = abs(i_normal * (centre_j - centre_i).unit());
									Real cosj = abs(j_normal * (centre_j - centre_i).unit());
									Real R_ij = (centre_j - centre_i).norm();

									if (_all_element[i]->_elem->dim() == 2)
										VF_ij += cosi * cosj * aera_j * aera_i / (pi * R_ij * R_ij);
									else if (_all_element[i]->_elem->dim() == 1)
										VF_ij += cosi * cosj * aera_j * aera_i / (2.0 * R_ij);
								}
							}
						}
						
						VF_ij = VF_ij / Ai;
					}
					_all_element[i]->_local_RD[j] = VF_ij;
				}
			}

			_communicator.barrier();
			for (unsigned int i = 0; i < all_element_size; i++)
			{
				_communicator.sum<Real>(_all_element[i]->_local_RD);
			}
			_communicator.barrier();

			cout << "Computing radiation view factor compeleted" << endl;

			//cout << endl;
			//cout << "/////////////////////Result of VF//////////////////////" << endl;
			//for (unsigned int ii = 0; ii < _all_element.size(); ii++)
			//{
			//	for (unsigned int i = 0; i < _all_element.size(); i++)
			//	{
			//		cout << _all_element[ii]->_local_RD[i] << " ";
			//	}
			//	cout << endl;
			//}

			RD_normalized(_is_seal);//归一化
			
			vector< vector<Real> > VF;
			vector< vector<Real> > RD_o;
			vector< vector<Real> > RD_n;
			vector<Real> eps_i;
			vector<Real> rho_i;

			eps_i.resize(all_element_size, 0.0);
			rho_i.resize(all_element_size, 0.0);
			for (int i = 0; i < all_element_size; ++i)
			{
				vector<Real> rd_temp;
				rd_temp.resize(all_element_size, 0.0);
				VF.push_back(rd_temp);
				RD_o.push_back(rd_temp);
				RD_n.push_back(rd_temp);
			}

			for (int i = 0; i < all_element_size; ++i)
			{
				eps_i[i] = _all_element[i]->_absorptivity;
				rho_i[i] = 1.0 - _all_element[i]->_transmissivity - eps_i[i];
				for (int j = 0; j < all_element_size; ++j)
				{
					VF[i][j] = _all_element[i]->_local_RD[j];
					RD_o[i][j] = _all_element[i]->_local_RD[j];
				}
			}

			for (int i_out = 0; i_out < _max_reflect_count; ++i_out)
			{
				cout << "Computing reflection [" << i_out + 1 << "]/[" << _max_reflect_count << "]" << endl;

				for (int i = 0; i < all_element_size; ++i)
				{
					for (int j = 0; j < all_element_size; ++j)
					{
						RD_n[i][j] = 0.0;
					}
				}

				_communicator.barrier();
				for (int i = local_start_id; i < local_end_id; ++i)
				{
					for (int j = 0; j < all_element_size; ++j)
					{
						RD_n[i][j] = eps_i[j] * VF[i][j];
						for (int k = 0; k < all_element_size; ++k)
						{
							RD_n[i][j] += VF[i][k] * rho_i[k] * RD_o[k][j];
						}
					}
				}
				_communicator.barrier();
				for (unsigned int i = 0; i < all_element_size; i++)
				{
					_communicator.sum<Real>(RD_n[i]);
				}
				_communicator.barrier();

				for (int i = 0; i < all_element_size; ++i)
				{
					for (int j = 0; j < all_element_size; ++j)
					{
						RD_o[i][j] = RD_n[i][j];
					}
				}
			}

			for (int i = 0; i < all_element_size; ++i)
			{
				for (int j = 0; j < all_element_size; ++j)
				{
					_all_element[i]->_local_RD[j] = RD_n[i][j];
				}
			}

			cout << "Computing radiation transfer factor compeleted" << endl;

			//////////////////////结束计算////////////////////////////
			_communicator.barrier();
		}

		else
		{
			ompComputeRD(_n_threads, _max_reflect_count);
		}
		
		Deal_with_BSMC(_BSMC);

		cout << endl;
		cout << "Writing......" << endl;

		ofstream rd_write;
		rd_write.open(_filename.c_str(),ios::trunc);
		if(!rd_write.is_open())
			mooseError("Error opening file '" + _filename + "' fail.");

		rd_write << _all_element.size() << " " << _all_element.size() << endl;

		for (unsigned int ii=0;ii<_all_element.size();ii++)
		{
			for (unsigned int i=0;i<_all_element.size();i++)
			{
				rd_write << _all_element[ii]->_local_RD[i] << " " ;
			}
			rd_write << endl;
		}

		rd_write.close();

		cout << "Write completed" << endl;
	}

	else
	{
		cout << endl;
		cout << "Reading......" << endl;

		_communicator.barrier();

		for (unsigned int ii=0;ii<_all_element.size();ii++)
		{
			_all_element[ii]->_local_RD.resize(_all_element.size(), 0);
		}

		ifstream rd_read(_filename.c_str());
		if(!rd_read.good())
			mooseError("Error opening file '" + _filename + "' from data.");

		string line;
		getline(rd_read, line);
		istringstream iss(line);
		int f;
		vector<int> head;
		while(iss >> f)
			head.push_back(f);

		unsigned int _num_i_of_RDij = head[0];
		unsigned int _num_j_of_RDij = head[1];

		if ( (_num_i_of_RDij != _all_element.size()) || (_num_j_of_RDij != _all_element.size()) )
			mooseError("Error file '" + _filename + "' connot march.");

		for (unsigned int i = 0; i < _all_element.size(); i++)
		{
			getline(rd_read, line);
			int j=0;
			istringstream transline(line);
			Real d;
			while (transline >> d)
			{
				_all_element[i]->_local_RD[j] = d;
				j+=1;
			}
		}

		_communicator.barrier();

		for(unsigned int i  = 0; i < _all_element.size(); i++)
		{
			_communicator.broadcast<Real>(_all_element[i]->_local_RD);
		}
		_communicator.barrier();

		cout << "Read completed" << endl;
	}

	end_time = clock();
	cout << endl;
	cout << "time: " << end_time-start_time << endl;

	//cout << endl;
	//cout << "/////////////////////Result of RD//////////////////////" << endl;
	//for (unsigned int ii = 0; ii < _all_element.size(); ii++)
	//{
	//	for (unsigned int i = 0; i < _all_element.size(); i++)
	//	{
	//		cout << _all_element[ii]->_local_RD[i] << " ";
	//	}
	//	cout << endl;
	//}
}

void ENRadiationHeatFluxUserObject::initialize()
{
}

void ENRadiationHeatFluxUserObject::execute()
{
	int findi = Find_i(_current_side_elem);
	Real temp_pow4_bar(0);
	Real temp_pow3_bar(0);
	Real temp_bar(0);

	for(unsigned int _qp = 0; _qp < _q_point.size(); ++_qp)
	{
		temp_pow4_bar += (_JxW[_qp]*pow(_temperature[_qp],4) );
		temp_pow3_bar += (_JxW[_qp]*pow(_temperature[_qp],3) );
		temp_bar += (_JxW[_qp] * _temperature[_qp]);
	}

	temp_pow4_bar /= _current_side_volume;
	temp_pow3_bar /= _current_side_volume;
	temp_bar /= _current_side_volume;
	_temperature_pow4_bar[findi] = temp_pow4_bar;
	_temperature_pow3_bar[findi] = temp_pow3_bar;
	_temperature_bar[findi] = temp_bar;
}

void ENRadiationHeatFluxUserObject::finalize()
{
	{
		_communicator.sum<Real>(_temperature_pow4_bar);
		_communicator.sum<Real>(_temperature_pow3_bar);
		_communicator.sum<Real>(_temperature_bar);
	}

	for(int i = 0; i < _all_element.size(); i++)
	{
		_flux_radiation[i] = 0;
		_flux_radiation_jacobi[i] = 0;
	}

	computeRadiationFlux();

	{
		_communicator.sum<Real>(_flux_radiation);
		_communicator.sum<Real>(_flux_radiation_jacobi);
	}

	for(int i = 0; i < _all_element.size(); i++)
	{
		_temperature_pow4_bar[i] = 0;
		_temperature_pow3_bar[i] = 0;
		_temperature_bar[i] = 0;
	}
}

void ENRadiationHeatFluxUserObject::ompComputeRD(int n_threads, int max_reflect_count)
{
	int all_element_size = _all_element.size();

	cout << "Computing radiation view factor start" << endl;
#pragma omp parallel for num_threads(n_threads) \
	schedule(static,1)
	for (int i = 0; i < all_element_size; ++i)
	{
		Real Ai = _all_element[i]->_elem->volume();
		int i_sep = _all_element[i]->_Disintegrate_centre.size();
		Point i_normal = _all_element[i]->getSideElementNormal();
		_all_element[i]->_local_RD.resize(all_element_size, 0.0);
		for (int j = 0; j < all_element_size; ++j)
		{
			Real VF_ij = 0.0;
			if (i == j)
			{
				VF_ij = 0.0;
			}
			else
			{
				Point j_normal = _all_element[j]->getSideElementNormal();
				int j_sep = _all_element[j]->_Disintegrate_centre.size();
				for (int i_in = 0; i_in < i_sep; ++i_in)
				{
					Point centre_i = _all_element[i]->_Disintegrate_centre[i_in];
					Real aera_i = _all_element[i]->_Disintegrate_aera[i_in];
					for (int j_in = 0; j_in < j_sep; ++j_in)
					{
						Point centre_j = _all_element[j]->_Disintegrate_centre[j_in];
						Real aera_j = _all_element[j]->_Disintegrate_aera[j_in];
						if ((_all_element[j]->getSideElementNormal()) * ((centre_j - centre_i).unit()) > (-RaylineZero))
						{
							continue;
						}
						else if ((_all_element[i]->getSideElementNormal()) * ((centre_i - centre_j).unit()) > (-RaylineZero))
						{
							continue;
						}
						else if (Is_Shelter_Engineer(i, j, centre_i, centre_j))
						{
							continue;
						}
						else
						{
							Real cosi = abs(i_normal * (centre_j - centre_i).unit());
							Real cosj = abs(j_normal * (centre_j - centre_i).unit());
							Real R_ij = (centre_j - centre_i).norm();

							if (_all_element[i]->_elem->dim() == 2)
								VF_ij += cosi * cosj * aera_j * aera_i / (pi * R_ij * R_ij);
							else if (_all_element[i]->_elem->dim() == 1)
								VF_ij += cosi * cosj * aera_j * aera_i / (2.0 * R_ij);
						}
					}
				}

				VF_ij = VF_ij / Ai;
			}
			_all_element[i]->_local_RD[j] = VF_ij;
		}
		//cout << i << " of " << all_element_size << endl;
	}
	cout << "Computing radiation view factor compeleted" << endl;

	RD_normalized(_is_seal);//归一化

	vector< vector<Real> > VF;
	vector< vector<Real> > RD_o;
	vector< vector<Real> > RD_n;
	vector<Real> eps_i;
	vector<Real> rho_i;

	for (int i = 0; i < all_element_size; ++i)
	{
		eps_i.resize(all_element_size, 0.0);
		rho_i.resize(all_element_size, 0.0);
		vector<Real> rd_temp;
		rd_temp.resize(all_element_size, 0.0);
		VF.push_back(rd_temp);
		RD_o.push_back(rd_temp);
		RD_n.push_back(rd_temp);
	}

	for (int i = 0; i < all_element_size; ++i)
	{
		eps_i[i] = _all_element[i]->_absorptivity;
		rho_i[i] = 1.0 - _all_element[i]->_transmissivity - eps_i[i];
		for (int j = 0; j < all_element_size; ++j)
		{
			VF[i][j] = _all_element[i]->_local_RD[j];
			RD_o[i][j] = _all_element[i]->_local_RD[j];
		}
	}

	for (int i_out = 0; i_out < max_reflect_count; ++i_out)
	{
		cout << "Computing reflection [" << i_out + 1 << "]/[" << max_reflect_count << "]" << endl;

		for (int i = 0; i < all_element_size; ++i)
		{
			for (int j = 0; j < all_element_size; ++j)
			{
				RD_n[i][j] = 0.0;
			}
		}
#pragma omp parallel for num_threads(n_threads) \
	schedule(static,1)
		for (int i = 0; i < all_element_size; ++i)
		{
			for (int j = 0; j < all_element_size; ++j)
			{
				RD_n[i][j] = eps_i[j] * VF[i][j];
				for (int k = 0; k < all_element_size; ++k)
				{
					RD_n[i][j] += VF[i][k] * rho_i[k] * RD_o[k][j];
				}
			}
		}

		for (int i = 0; i < all_element_size; ++i)
		{
			for (int j = 0; j < all_element_size; ++j)
			{
				RD_o[i][j] = RD_n[i][j];
			}
		}
	}

	for (int i = 0; i < all_element_size; ++i)
	{
		for (int j = 0; j < all_element_size; ++j)
		{
			_all_element[i]->_local_RD[j] = RD_n[i][j];
		}
	}
	cout << "Computing radiation transfer factor compeleted" << endl;
}

void ENRadiationHeatFluxUserObject::Deal_with_BSMC(int BSMC)
{
	for(int k = 0; k < BSMC; ++k)
	{
		vector<Real> N_i;
		N_i.resize(_all_element.size(), 1.0);
		Real temp = 0;
		for (unsigned int i=0;i<_all_element.size();i++)
		{
			Real normalized_parameter = _all_element[i]->_elem->volume()*_all_element[i]->_absorptivity;
			N_i[i] = N_i[i]*normalized_parameter;
			for (unsigned int j=0;j<_all_element.size();j++)
			{
				_all_element[i]->_local_RD[j] = _all_element[i]->_local_RD[j]*normalized_parameter;
			}
		}
		for (unsigned int i=0;i<_all_element.size();i++)
		{
			for (unsigned int j=0;j<_all_element.size();j++)
			{
				N_i[i] += _all_element[j]->_local_RD[i];
			}
		}
		for (unsigned int i=0;i<_all_element.size();i++)
		{
			temp = 0;
			for (unsigned int j=i;j<_all_element.size();j++)
			{
				temp = _all_element[i]->_local_RD[j] + _all_element[j]->_local_RD[i];
				_all_element[i]->_local_RD[j] = temp/N_i[i];
				_all_element[j]->_local_RD[i] = temp/N_i[j];
			}
		}
	}
}

void ENRadiationHeatFluxUserObject::computeRadiationFlux()
{
	Real sigma = 5.67e-8;
	Real Flux_Rad = 0.0;

	int local_begin, local_end, local_n;
	int quotient = _all_element.size()/n_processors();
	int remainder = _all_element.size()%n_processors();
	if(processor_id() < remainder)
	{
		local_n = quotient + 1;
		local_begin = processor_id()*local_n;
		local_end = local_begin + local_n;
	}
	else
	{
		local_n = quotient;
		local_begin = processor_id()*local_n + remainder;
		local_end = local_begin + local_n;
	}

	for (int i = local_begin; i < local_end; i++)
	{
		Flux_Rad = 0.0;
		for (unsigned int j = 0;j < _all_element.size(); j++)
		{
			Flux_Rad += (_all_element[j]->_local_RD[i])*(_all_element[j]->_elem->volume())*_all_element[j]->_absorptivity*_temperature_pow4_bar[j];
		}
		_flux_radiation[i] = sigma*Flux_Rad/_all_element[i]->_elem->volume();
		_flux_radiation_jacobi[i] = 4*sigma*(_all_element[i]->_local_RD[i])*_all_element[i]->_absorptivity*_temperature_pow3_bar[i];
	}
}

int ENRadiationHeatFluxUserObject::Find_i(const Elem * elem) const
{
	int findi = -1;

	for (unsigned int i = 0; i < _all_element.size(); i++)
	{
		if( (_all_element[i]->_elem->centroid() - elem->centroid()).norm() < RaylineZero)
		{
			findi = i;
			break;
		}
	}

	if(findi != -1)
		return findi;
	mooseError("Cannot find i.");
}

bool ENRadiationHeatFluxUserObject::Is_Shelter_Engineer(int i_index, int j_index, Point i_centre, Point j_centre)
{
	int sideelement_size = _all_element.size();
	int j_wanted = -1;
	RayLine ray = RayLine(i_centre, j_centre, 1);
	Point point_temp = ray._p0;
	Point p1 = ray._p1;

	for (int j = 0; j < sideelement_size; ++j)
	{
		if (j == i_index || j == j_index)
		{
			continue;
		}
		else if ((_all_element[j]->getSideElementNormal()) * (ray._normal) > (-RaylineZero))
		{
			continue;
		}
		else if (!(ray.sideIntersectedByLine(_all_element[j]->_elem, point_temp)))
		{
			continue;
		}
		else if ((point_temp - ray._p0).norm() < (p1 - ray._p0).norm())
		{
			return true;
		}
		else
		{
			continue;
		}
	}
	return false;
}

void ENRadiationHeatFluxUserObject::RD_normalized(bool is_seal)
{
	if (is_seal)
	{
		int all_element_size = _all_element.size();
		for (int i = 0; i < all_element_size; ++i)
		{
			Real Summation = 0.0;
			for (int j = 0; j < all_element_size; ++j)
			{
				Summation += _all_element[i]->_local_RD[j];
			}
			for (int j = 0; j < all_element_size; ++j)
			{
				_all_element[i]->_local_RD[j] = _all_element[i]->_local_RD[j] / Summation;
			}
		}
	}
}

void ENRadiationHeatFluxUserObject::print_mesh_to_file()
{
	string _filename = "mesh_for_fast_flux.dat";

	ofstream rd_write;
	rd_write.open(_filename.c_str(), ios::trunc);
	if (!rd_write.is_open())
		mooseError("Error opening file '" + _filename + "' fail.");

	rd_write << _all_element[0]->_elem->dim() << '	' << _all_element.size() << endl;
	for (unsigned int ii = 0; ii < _all_element.size(); ii++)
	{
		Point p = Point(0.0, 0.0, 0.0);
		int n_nodes = 0;
		if (_all_element[ii]->_elem->type() == TRI3)
		{
			n_nodes = 3;
		}
		else if (_all_element[ii]->_elem->type() == QUAD4)
		{
			n_nodes = 4;
		}
		else
		{
			mooseError("Error! This sideelement is not TRI3 or QUAD4. Please check it!");
		}
		rd_write << n_nodes << endl;
		for (unsigned int jj = 0; jj < n_nodes; jj++)
		{
			p = _all_element[ii]->_elem->point(jj);
			rd_write << p(0) << '	' << p(1) << '	' << p(2) << endl;
		}
		p = _all_element[ii]->getSideElementNormal();
		rd_write << p(0) << '	' << p(1) << '	' << p(2) << endl;
		rd_write << _all_element[ii]->_transmissivity << '	' << _all_element[ii]->_absorptivity << endl;
		rd_write << _all_element[ii]->_transmissivity << '	' << _all_element[ii]->_absorptivity << '	'
			<< 1.0 << '	' << 0.0 << endl;
	}

	rd_write.close();

	cout << "Write elem_T_list to " << _filename << " completed" << endl;
}

void ENRadiationHeatFluxUserObject::print_T_to_file()
{
	string _filename = "T_for_fast_flux.dat";
	Point p = Point(0.0, 0.0, 0.0);

	ofstream rd_write;
	rd_write.open(_filename.c_str(), ios::trunc);
	if (!rd_write.is_open())
		mooseError("Error opening file '" + _filename + "' fail.");

	rd_write << _all_element.size() << endl;
	for (unsigned int i = 0; i < _all_element.size(); i++)
	{
		p = _all_element[i]->_elem->centroid();
		rd_write << p(0) << '	' << p(1) << '	' << p(2) << '	' << _temperature_bar[i] << endl;
	}

	rd_write.close();
}

void ENRadiationHeatFluxUserObject::print_flux_to_file()
{
	string _filename = "receive_flux.dat";
	Point p = Point(0.0, 0.0, 0.0);

	ofstream rd_write;
	rd_write.open(_filename.c_str(), ios::trunc);
	if (!rd_write.is_open())
		mooseError("Error opening file '" + _filename + "' fail.");

	rd_write << _all_element.size() << endl;
	for (unsigned int ii = 0; ii < _all_element.size(); ii++)
	{
		p = _all_element[ii]->_elem->centroid();
		rd_write << p(0) << '	' << p(1) << '	' << p(2) << '	' << _flux_radiation[ii] << endl;
	}

	rd_write.close();
}
